Coulomb gap in a model with finite charge transfer energy. 
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The Coulomb gap in a donor-acceptor model with finite charge transfer energy A describing the 
electronic system on the dielectric side of the metal-insulator transition is investigated by means of 
computer simulations on two- and three-dimensional finite samples with a random distribution of 
equal amounts of donor and acceptor sites. Rigorous relations reflecting the symmetry of the model 
presented with respect to the exchange of donors and acceptors are derived. In the immediate 
neighborhood of the Fermi energy the the density of one-electron excitations g(e) is determined 
solely by finite size effects and g(e) further away from fj, is described by an asymmetric power law 
with a non-universal exponent, depending on the parameter A. 



PACS numbers: 71.23.-k, 71.30. +h, 71.45. Gm 
I. INTRODUCTION 

Doping of solids might lead to drastic qualitative 
changes in their properties. The metal- insulator tran- 
sition (MIT) is a spectacular manifestation of this. The 
understanding of the driving forces of the MIT is a long=. 
standing problem. In the early seventies, the prediction!!! 
was made that on the dielectric side of the MIT the long- 
range Coulomb interactions deplete the density of one- 
electron excitations (DOE) g{e) near the Fermi energy 
\i. Further, analytical calculations of g{e) with Coulomb 
correlation taken into consideration have been performed 
on the metallic side of the MIT. Altshuler and AronovEl 
showed that for the metallic case g(e) in three dimensions 
has a cusp-like dependence g(e) ~ |e — ii] 1 ^ 2 near /i. This 
was later confirmed in electron tunneling experiments for 
amorphous alloys! and granular metalsQ. 

On the insulating side of the MIT charge transport 
occurs via inelastic electron tunneling hopping between 
states localized on the impurity sites with one-electron 
energies close to [i. Mottfl demonstrated that at low 
temperatures electrons seek accessible energy states by 
hopping distances beyond the localization length, lead- 
ing to a hopping conductivity u(T) ~ exp(— Tq/T) v with 
To being a characteristic temperature depending on lo- 
calization length and with the hopping exponent v = 1/4 
for the non-interacting case in three dimensions. Efros 
and ShklovskiiH (ES) argued that the ground state of a 
system with long-range Coulomb interactions is stable 
with respect to one-particle excitations only if g(s) in 
the vicinity of /i has the symmetric shape 

g^-le-^- 1 (1) 

with the universal exponent D — 1 depending only on the 
dimensionality D of the system. In particular, ES pre- 
dicted that in D = 3 g{e) — & (^) 3 (e — it) 2 , where \ 
is the dielectric constant and e is the electron charge. 



Because g[e) vanishes only at e = fj,, this is called 
a "soft" Coulomb correlation gap with a width Ae ~ 
e 3 (A /x 3 ) 1 ^ 2 , wherg N is the DOE far away from //. The 
power law ([!]) givesQ a hopping exponent v = D/(D + 3) 
at low temperatures, so for three-dimensional system 
with long-range Coulomb interactions v = 1/2. 

The intriguing hypothesis about universality of (]]]) has. 
stimulated furth&t^theoretical research, both analyticalQ 
and pimericalclt 3 ]. To establish the hypothesis ([j]) 
EfrosEfJ used the ground-state stability conditions for lo- 
calized electrons (LES) with respect to charge transfer 

e 2 

Ej -Ei > 0, (2) 

where £j and £j are the one-particle energies of a neutral 
donor on a site i and of a charged donor on a site j, re- 
spectively, and Tij is the distance between the sites i and 
j. The conditions (0) were used to heuristically derive 
a non-linear integral equation for and then as- 

symptotic analysis of this equation leadstfJ to the power 
law (0). 

LES have been studied using the so-called classical 
donor-acceptor {d-d) model (see, e. g. Ref. [l5|). Within 
this model, the system considered is modeled by a con- 
tinuous sample with randomly distributed k x N (k < 1) 
acceptor and N donor sites. Each acceptor site is neg- 
atively charged whereas out of N only k x N donors 
have a positive charge which leads to a large number 
of configurations of charged donors. Moreover, each of 
these configurations must obey not only conditions (|^) 
but also more complicated conditions related to many- 
particles excitations (e. g., charge transfer involving four, 
six, etc. sites). Efros conjecture about the universality 
implies that g(e) does not depend on peculiarities of the 
particular model aiuL-as-a consequence, further theoret- 
ical studies of LESlaOii 3 ] were confined to a lattice d-a 
model proposed in Ref. 113. In this model, N donors are 
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localized on all the sites of a D-dimensional lattice and 
the negative charge from k x N acceptors is uniformly 
smeared over the lattice sites so that each site i has a 
charge e(rij — fc), where ni = 1 if a donor on the site i is 
ionized and rii — if a donor is neutral. Disorder in this 
model is ensured by introducing randomly distributed 
one-site potentials. Monte Carlo simulationsc3 on very 
large specimens of the lattice d-a model, however, have 
given rise to doubts about the universality of the g(e) 
behavior. 

Another hint about possible non-universal behavior of 
g(e) has come from the intriguing and still not com- 
pletely unfolded problem whether the so called spin-glass 
phase does exist in the classical d-a model (see, e. g. Ref. 
|l8|-po|). Grannan and YuEa studied the classical three- 
dimensional d-a model with k = 0.5 but with the total 
acceptor charge uniformly distributed over donor sites 
as in the lattice d-a model. In this case, the classical d-a 
model is equivalent to a model of Ising spins, localized on 
randomly distributed sites, with pairwise Coulomb inter- 
actions, a modeLi:! which a transition into the spin-glass 
state was foundcS to occur at non-zero temperature. It 
was then concluded that such a transition should exist 
in all d-a models (with and without smearing of negative 
charge, defined on a lattice or on a continuous sample) as 
well because .of the Efros universality hypothesis. Voita 
and SchreberEa, however, have shown that the spin glass 
transition does not exist in the lattice d — a modeltj. 
Besides, in recent work by one of usa it was unequivo- 
cally demonstrated that the ground state of the classical 
d-a model and that of the model studied in Ref. [l8| are 
qualitatively different. An analysis of histograms H.[Q a p] 

of the so called overlaps Q a p = X^<K n ?' n f) (here a 
and (3 refer to different pseudo-ground states (PGS) ob- 
tained by direct descents) has revealed that, indeed, for 
the model studied in Ref. [l8| H[Q a f3] has a symmetric 
Gaussian shape with the maximum at {Qap) — and 
with the dispersion (Q a p) ~ N^ 1 . This means that a 
large number of microscopically different PGS's doss ex- 
ist in the model and according to Parisi's theorytil this 
implies the existence of a spin-glass state at low tem- 
peratures. Eurthcr Monte Carlo simulations at finite 
temperatureso revealed the typical finite-size scaling of 
the spin-glass susceptibility. In the classical d-a model, 
however, H[Q a p] has its maxima at {Q a p) — 1 which 
means that all PGS's generated are the same from mi- 
croscopical point of view. The absence of microscopically, 
different PGS's in the classical d-a model was explained^ 
by the pinning of all PGS's on the electric field created 
by the discretely distributed acceptor charges. 

Therefore, it is highly desirable to study the proper- 
ties of not only the classical d-a model, but of its various 
modifications as well. In the present work we consider a 
modified classical d-a model (MCDAM) in which accep- 
tors can be neutral, so the energy A of the charge transfer 
from a donor to an acceptor (d°+a° — > 



neutral acceptor, a charged acceptor, respectively) has to 
be finite. The classical d-a model might be then viewed 
as the limit of the MCDAM as A — > oo. We have inves- 
tigated the shape of the Coulomb gap (i. e. g(e) both for 
donor and acceptors in the vicinity of the Fermi level) in 
two- and three-dimensional MCDAMs at T=0 and found 
that the behavior of g(e) is in strong contradiction to 
the Efros conjecture about the universality of g(e). The 
rest of the paper is organized as follows. In Section II 
we introduce the MCDAM and arrive at some rigorous 
results which follow from a symmetry of the MCDAM 
with respect to the exchange of donor and acceptor sites. 
Further, the algorithm of energy minimization for the 
MCDAM including a discussion about inherent finite size 
effects is presented in Section III. Section IV is devoted 
to a description of the main results obtained. In Section 
V we discuss possible causes of universality violation in 
the MCDAM, analyze experimental data available in the 
literature and predict possible experimental situations in 
which the non-universal behavior of 17(e) might be ob- 
served. And finally, a summary is presented in Section 
VI. 



II. BACKGROUND 
A. Model 

We consider a D-dimensional system of volume L D , in 
which an equal number of acceptor and donor sites N 
are allocated according to the Poisson distribution with 
a density n — N x L~ D . It is convinient to choose en- 
ergy unit Eq as an energy of the Coulomb interaction 
between a pair of acceptors, say, localized on the average 
distance rT 1 ! , Eq = e 2 ?! 1 / 15 j\- In typical bulk semi- 
condustors n ~ 10 18 cm~ 3 and \ ~ 10; so Eo ~ 0.02 eV. 
Hereafter all expressions will be written in dimensionless 
units nT 1 / for length and Eq for energies. A micro- 
scopic state of a particular spatial arrangement of the 
donor and acceptor sites (henceforth referred to as the 
sample R) is determined by a set of occupation numbers 
(n a ,n d ) = {n a (i),n d (k), i = 1, 2, . . . , N, k = 1, 2, . . . , N} 
determined in the following way. For the acceptors, 
n a (i) = 1 if an acceptor on an acceptor site i has cap- 
tured an electron and n a (i) = if an acceptor is neutral. 
For the donors, n d (k) = 1 if a donor on a donor site k is 
neutral and n d {k) = if a donor has given an electron 
away. We investigate the LES from the dielectric side of 
the MIT, so ocb < 1 (&b is the localization length of the 
electron on donor). The energy of the sample, assuming 
that all the interactions are of Coulomb origin, then is 
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a , a~ stand for a neutral donor, a charged donor, a 
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where indices i,j and k,l number acceptor and donor 
sites, respectively, r^"° , and r^T are the distances 
between the acceptors on the sites i and j, between the 
donors on the sites k and I, and between the acceptor on 
the site i and the donor on the site k, correspondingly, 
and A is the energy of charge transfer between acceptor 
and donors. When charge transfer occurs in the system, 
the energy of the sample changes by 



SE(n a ,n d ) = ^e a {i)Sn a (i) + } j e d (k)Sn d (k) 
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where e a (i) is the one-electron excitation (OEE) energy 
for the acceptors 
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e f /(i) is the corresponding OEE energy for the donors 
and Sn a (i) (Sn d (k)) denotes the change of the occupa- 
tion number on the acceptor (donor) site i (fc). If a mi- 
croscopic state (n°,n°) of the sample is the ground-state 
of this sample then for any excitation the relation 



SE(n° a ,n° d ) >0 



(6) 



holds. The specific appearance of the conditions (||) de- 
pends on what excitations are allowed in the model sys- 
tem considered. 

In the present paper we investigate the simplest case 
when only pairs of sites are involved in the charge trans- 
fer which, in turn, is allowed to occur in four different 
ways: (i) via electron hops between a pair of the ac- 
ceptors {n a (i) = l,n a (j) = 0} -> {n a (i) = 0,n a (j) = 
1}; (ii) via electron hops between a pair of donors 
{n d (k) = l,n d (l) = 0} — {n d (k) = 0,n d (l) = 1}; 
(iii) via ionization process {n a (i) = 0,n d (k) — 1} — > 
{n a (i) = l,n d (k) — 0} and (iv) via recombination pro- 
cess {n a (i) = l,n d (k) = 0} -> {n a (i) = 0,n d (k) = 1}. 
For each of those processes there is an unique set of 
{8n a (i), Sn a (j), 5n d (k), Sn d (l)}. For instance, for the 
acceptor- acceptor hops 

Sn a (i) = -1, 5n a (j) = 1, 5n d (k) = 0, 8n d {l) = 0. (7) 

Substituting (|7|) into (||) one obtains the ground-state 
stability relation with respect to the charge transfer be- 
tween the pair of acceptors on the sites i and j 



where e^°\i) denotes e a (i) if n(i) — 1(0). The stability 
conditions with respect to the other three manners of the 
charge transfer are obtainable in the similar manner. 

The relation (||) implies that e Q 's for the neutral ac- 
ceptors are, in general, larger than e a 's for the charged 
acceptors. Furthermore, the pair of neutral and charged 
acceptors might be located on any distance and therefore 
in the thermodynamic limit the chemical potential for the 
acceptors (i. e. an energy level which separates the ener- 
gies of the neutral and charged acceptors) is determined 
as 
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i{ E °(i)}=max{4(i)}. 



(9) 



Alike, there exist the chemical potential fi d for the donors 
as well. Moreover, the stability relations with respect to 
the ionization and recombination lead to 



/ia = Md = M- 



(10) 



Despite the finite size of samples we investigated, the re- 
lation ([n]) with the chemical potentials calculated from 
@ is valid within the limits of accuracy of our calcula- 
tions (see Sect. III). 

A macroscopic state of the sample R is characterized 
by degree of acceptor ionization 



by the DOE for acceptors 
1 ■ 



g a (e a ,R) 



N ^ 



5(e - e a (i)) 



(11) 
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and by the corresponding DOE g d (e d ,Il) for the donors. 
Note, that for the finite samples (especially for the rel- 
ative small systems we were able to investigate) C a (R), 
5a(e a ,R) and g d (e d ,H) depend essentially on the partic- 
ular implementation R of the spatial distributions of the 
donor and acceptor sites (if a sample would be big enough 
all quantities would be self-averaging). Therefore, in 
order to obtain reliable results, one has to work with 
the quantities C a = (C Q (R)), g a (e) = (p (e ,R)) and 
g d (e) = (g d {e d , R)), where (. . .) denotes the average over 
a number of R's. Note, that the values <7a(d)( £ a(d)i R)<fe 
obtained for independent R's are scattered according to 
the Gaussian distribution with the mean g a (d){ £ )de and 
the standard deviation ^ g a ( d ) (e)de. In the region of the 
Coulomb gap g a ( d ){e)d£ ~ 10~ 4 and dispersion is several 
orders of magnitude larger than the mean. Therefore, in 
order to reduce the statistical noise in the final g a (d){ £ ) 
dependences an average is needed over a sufficient large 
amount of independent samples (we performed calcula- 
tions with up to 10 samples). 
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B. Acceptor-donor symmetry 



III. METHOD 



Let us rewrite the energy (|||) in terms of the OEE 
energies (|j) 

E(n a ,n d ) = - ^2e a (i)n a (i) - 

i 

-jE^ki-^-jE^. (13) 

k i 

The system investigated is electrically neutral, i. e. for 
any sample 
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(l-n d (k)). 



(14) 



Then, the energies of the microscopic states (n , n d ) and 
(n a ,n d ) for a sample R and its "mirror" reflection R* 
(when the donor and acceptor sites exchange places keep- 
ing the spatial arrangement of sites unchanged) , are equal 
under the following conditions 



e a (i)+e* d (i) = s d (k) + s* a (k) 



-A 



(15) 



and 



<(*) = 0--na(i)) 



i* d (k) = (l-n a (k)). (16) 



The stability relations (@) for the ground-state (n a ,n d ) 
of the sample R transform into stability relations for the 
ground-state (n a *,n d *) of the sample R* through the re- 
lations (|l5|,|l6]) as well. 

Since averaging over samples includes all possible pairs 
R and R*, it follows from the symmetry relations Jl5| ) 
and ( |l6| ) along with the definition (12) that g a (e) can be 
mapped to g d {e) using the relation 



9d{e) =g a (-e-A) 



(17) 



The symmetry of the model imposes also a relation 
between the Fermi energy /i (|9pC|) and the parameter 
A of the model. Expressing n a ^(i[k]) in terms of the 
Heaviside's step functions n a u}\ = #(/i — e a [ d ](i[k])), the 
quantity C a can be written in the form 



C a 



g a (e)de 



9d{e)de 



(18) 



The symmetry relation ( |17| ) transforms ( |18|) into an in- 
tegral relation 



g d {e)de = / g d (e)de, 

— — A J 



which has a meaning only if 

A 



(19) 



(20) 



Thus, the Fermi energy of our model system in the ther- 
modynamic limit is a fundamental quantity depending 
only on the energy of charge transfer from an acceptor 
to a donor. 



A. Algorithm of energy minimization 

We start from a random allocation of N donor and 
N acceptor sites in the continuous D-dimensional sys- 
tem (generate a sample R) with the density n = 1, so 
that the system has a linear size L = N 1 / and then 
charge randomly chosen C a x N both donors and accep- 
tors (usually we take C a = 0.7), i. e. generate an initial 
microscopic state (IMS) (n a ,n d ) of the sample R. Fur- 
ther, we search for such microscopic state (n a ,n d ) which 
obeys the stability conditions (||) with respect to the four 
mechanisms of the charge transfer allowed in our model. 
We used an algorithm which is an extension of the al- 
gorithm proposed in Ref. || to the case A ^ oo. The 
algorithm consists of the three main steps. 

In order to save computer time, first, we look for pairs 
a — a~ (d° — d + ) for which the "crude" stability relation 
Ae 



£ 







a(d) c a(d) 



> is violated. Then, the energy of 
the system is decreased by transferring an electron be- 
tween such pair of sites for which Ae has its minimal 
non-positive value. This process is repeated until a state 
is reached, in which Ae > for all possible a — a~ and 
d° — d + pairs (step I). In the similar manner, we further 
minimize the energy of the system with respect to the 
"true" stability relations (||) for the charge transfer be- 
tween the a° — a~ and d° — d + pairs (step II). And, finally, 
in the step III we diminish the energy of the system with 
respect to the stability relations for ionization and recom- 
bination processes. Since ionization and recombination 
processes change the degree C a of the system ionization, 
each time after one of these processes takes place dur- 
ing calculations, we go back to the step II. Repeating 
the steps II and III, we finally arrive at a microscopic 
state (n a ,n®) for which all four stability conditions are 
fulfilled. We name the procedure (n a ,n d )^ {n a ,n° d ) via 
above steps I, II and III as "a single descent" . 

It should be noted, however, that the state (n a ,n d ) is 
not necessarily the ground state of the sample R since 
for the ground state, in general, not only the simplest re- 
lations (|8j) with only pairs of sites included, but the more 
complicated relations involving quadruplets, sextets, etc. 
of sites have to be fulfilled. Therefore, the state (n a ,n d ) 
(after Ref. |9|) hereafter will be referred to as the pseudo- 
ground state (PGS) of the sample R. Then, two questions 
naturally arise: How close the PGS and the ground state 
of the given sample are and how this may influence the 
output of our calculations? In order to answer the first 
question, we calculate and analyze the histograms TL for 
the so-called overlaps 



(21) 



where indices a and (3 refer to PGS's which are obtained 
by means of the single descent on the same sample but 
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with different IMS (n a ,nd). If two PGS's are identical 
then Q aj 3 = 1. We calculated for the D = 2 system 
with N = 500 at A = the mean Q(R) = {Q af )) a p for 
the sequence of 100 PGS's generated by single descents 
from the different IMS of the same sample R. We further 
acquire Q(R) for 100 different samples and obtain that 
the mean Q = (Q(R))r = 0.96. It means that in PGS 
generated by the single descent only 20 acceptors out of 
500 are, in average, in the "wrong" states compared to 
those in the true ground state of the sample. 

In order to evaluate how the "erroneousness" of PGS 
influences the outcome of our calculations we perform an 
analysis of ground states obtained by means of the so 
called multirank descents. Descent of rank to comprises 
of a consequence of the single descents on the same sam- 
ple with different IMS when calculations are stopped af- 
ter the lowest observed PGS energy repeats to times. We 
calculate Q (all other parameters were the same as de- 
scribed in the previous paragraph, where actually the 
case to = was explored) for descents with different 
ranks to = 5, 10, 15 and found that, for instance, for 
to = 15 (which implies drastic increase in the compu- 
tation time) Q = 0.990. g a (s) and gd(s) obtained from 
the PGS's generated by means of the single descents and 
by means of descents with to = 10, say, do not differ 
within the limits of statistical errors. So, we conclude, 
that reliable results can be obtained by means of single 
descents already, thereby saving a lot of computer time 
and resources. 



B. Finite-size effects 

Due to constraints in computer resources, the largest 
samples, we were able to deal with, comprise up to N — 
2000 donor and N = 2000 acceptor sites (L ~ 45 for 
D = 2 and L ~ 12 for D = 3). Such relative small sizes 
of the samples investigated might influence the outcome 
of calculations. Detailed analysis of finite size effects on 
the results obtained will be presented in Section IV and 
here we want to make two remarks about inherent finite 
size effects in the model system considered. 

First, as follows from (H), the energies e° for the neutral 
acceptors and e\ for the charged ones in finite samples 
at T = cannot be further away than (L x \/~D)~ 1 . This 
implies that g{e a ) = within the e a interval 

\e a - fj,\ < (2L x VD)- 1 (22) 

Of course, the same holds for donors as well. The rela- 
tion ( |22| ) gives the estimation how close to /i data on the 
energy spectrum are, in principle, obtainable from the 
calculations on finite samples. 

Secondly, as follows from (||) the energies e a and Ed 
for the finite samples are sensitive to the location of the 
donor and acceptor sites. Therefore, the Fermi energy fj, 
for finite samples does differ, in general, from sample to 



gfe-n) 

0.0004 



0.0002 





-0 

FIG. 1. Density of one electron excitations g a (e — fj,) 
in the vicinity of the Fermi energy fj, obtained for the 
two-dimensional model with N = 1500 at A = (cir- 
cles), 2 (squares), 4 (diamonds) and 10 (triangles). Data 
points presented in the figure are calculated as the average 
over 10.000 (A = 0), 5.100 (A = 2), 3.700 (A = 4) and 2.200 
(A = 10) different samples. Insert shows double logarithmic 
plot of g a (e — /i) for e > fj, in the region e — fi < 0.05. 

sample. A straightforward averaging of g(e) over differ- 
ent samples might thus lead to a distortion of the g(e) 
shape especially in the region where the Coulomb gap is 
observed. In order to avoid this undesired effect, we used 
a trick first proposed in Ref. ^|. During accumulation of 
the results for g(e) we added together g(s) for the same 
values of e — /x(R) rather than for the same values of e. 
Here /x(R) denotes the Fermi energy for a finite sample 
R calculated as 

M (R) = \ (min^WI + max^)}), (23) 

Such way of doing g{e) average entirely excludes the in- 
fluence of the fluctuations of the Fermi energy in the finite 
samples on the shape of the Coulomb gap. 

Finally, we remark that all the data presented below 
were obtained for the open boundary condition. In order 
to ensure that results obtained are not determined by the 
type of the boundary conditions used in calculations, we 
performed calculations of the two-dimensional MCDAM 
at A = with different N and found that periodic bound- 
ary conditions only effectively reduce the linear size of a 
sample, leaving the qualitative shape of the parameters 
calculated unchanged. 

IV. RESULTS 

According to the symmetry relation ( |l7| ) g a (s) and 
<7d(e) can be easily mapped to each other for any val- 
ues of e and hence all the results presented below con- 
cern the acceptor sites only. One can expect that the 
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FIG. 2. Density of one electron excitations g a (s — /i) for 
£ > jJ. (a,b) and e < [i (c,d) obtained for the two-dimensional 
model @ at A = (a,c) and 4 (b,d), with N = 500 (curves 
numbered 1), 1000 (2) and 1500 (3). The dashed lines are 
least-squares power-law fits g a (e — /x) ~ |e — /i| 7 with 7 = 0.9 
(a), 0.55 (b), 0.98 (c) and 0.78. Data presented in the figure 
are calculated as the average over 10.000 different samples 
(except the case N = 1500 and A = 4 with the average over 
3700 different samples). 

width of the Coulomb gap Ae and the energy scale in 
our model Eq = e 2 ™}/ j\ are of the same order of mag- 
nitude. Fig.[l] shows g a (e — /j) in the vicinity of the Fermi 
energy /1 obtained for the two-dimensional samples with 
N = 1000 and various values of A. As it is seen, g a {e — n) 
depends considerably on A except for a narrow window 
\e — fx\ < 0.05, where all data merge into some "uni- 
versal" curve symmetric with respect to the curve 
which can be anticipated to obey the Efros universal- 
ity hypothesis (|l|). However, a double- logarithmic plot 
of the "universal" g a (s — fi) (insert in the Fig .[[J), reveals 
that the behavior of g a (s — fi) in the "universality" region 
is not even a power law. The width of this "universality" 
region is comparable to the width of the region where 
g a (e — /i) = due to the finite size effects (for the data 
presented in Fig.|l| relation ( 122 ) gives |e — jtij < 0.011), so 
it is plausible to suggest that the "universal" behavior of 
g a (e — /i) is governed by the finite-size effects. This is 
clearly demonstrated in Fig.|| where g a (e — pb) are shown 
for several sizes of the samples investigated. 

The e window where finite size effects are severe, 
shrinks considerably with increasing N for all values of 
A we investigated. For instance, g a (e — fi) for N = 500 
and N = 1000 at A = (see Fig.^a,c) merge when 
\e— fi\ > 0.2 while corresponding curves for N = 1000 and 
N = 1500 arc indistinguishable already at \e — fi\ > 0.1. 
The statistical noise observed for the curves in Figj^ is 
quite small even close to fi and hence, the influence of 
insufficient large statistics on the results obtained is ex- 
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FIG. 3. The exponent 7 of the power law g a (e— A 4 ) ~ l e ~ 
as a function of the charge-transfer energy A. The data 
are obtained from least-squares fits of g a (£ — A 4 ) for the 
two-dimensional model ([]) with iV = 1500 within the region 
0.2 < |e — fi\ < 0.7. Circles represent the positive values of 
e — n while diamonds stand for the negative values of e — fi. 
Lines are guides to the eye. 

eluded. Note, that the "universal" behavior of g(e) in 
the vicinity of [i obtained for the classical d — a model 
(see Fig. 3 in Ref. |ll|) is most likely due to the finite size 
effects as well. 

In the region |e — /i| > 0.2, where the curves for all 
N collapse into a single curve (and where we believe the 
thermodynamic limit is reached), the behavior of g a (s — 
/i) is described by a power law g a (e — /1) ~ \e — /i| 7 . The 
deviation from the power-law observed far away from fi 
(|e — > 0.7) is due to the boundaries of the Coulomb 
gap which, as was mentioned above, are ~ 1 in units of 
Eq. One can see from a comparison of the data shown 
in Fig.|| for different A, that the exponent 7 depends 
considerably on A. Furthermore, values of 7 in the region 
e — [i > and those in the region e — [i < differ as well 
with this difference increasing with increasing A. The 
data for 7 obtained for the two-dimensional MCDAM 
are summarized in Fig.|| where a significant deviation of 
7 from the value D — 1 predicted by the hypothesis (|l|) 
is observed at all values of A investigated except for the 
case A = when 7 w 1 within the limits of statistical 
accuracy. Note, that the deviation of 7 from its predicted 
value grows monotonically with increasing A. At A = 10 
where the features of the MCDAM are expected to be 
nearly the same as those of the classical d — a model 
with all the acceptors being ionized (indeed, the degree of 
the acceptor ionization C a ~ 0.9 for the two-dimensional 
MCDAM at A = 10, see Fig.| below) the deviation from 
the Efros exponent is very large. 

The main results for g a (e — ji) obtained for the three- 
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FIG. 4. Density of one electron excitations g a (e — ft) 
in the vicinity of the Fermi energy /i obtained for the 
three-dimensional model (^) with N = 1000 at A = (cir- 
cles), 2 (squares), 4 (diamonds) and 10 (triangles). Data 
points presented in the figure are calculated as the av- 
erage over 10.000 different samples. Inserts show dou- 
ble-logarithmic plots of g a (e — fj.) at A = 2, for N = 500 
(curves numbered 1), 1000 (2) and 2000 (3), in the regions 
e > fi (a) and for e < \jl (b). The dashed lines in the in- 
serts are least-squares power-law fits g a (s — \i) ~ e — /i| 7 with 
7 = 1.16 (a), 1.29 (b), 




FIG. 5. The exponent 7 of the power law 

g a {E — [h) ~ |e — /i| 7 as a function of the charge-transfer energy 
A. The data are obtained from least-squares fits of g a (s — fi) 
for the three-dimensional model with N = 1000 within 
the region 0.4 < \s — fi\ < 0.8. Circles represent the positive 
values of s — \i while diamonds stand for the negative values 
of e — fi. Lines are guides to the eye. 




FIG. 6. The degree of acceptor ionization C a as a function 
of the charge-transfer energy A. The data are obtained for the 
model (^) in two (circles) and three (diamonds) dimensions 
with N = 500 as an average over 1000 different samples. The 
solid lines are third-degree polynomial fits. 



dimensional MCDAM are summarized in Figs. || and ||. 
It is seen, that the behavior of g a (e — fj.) in three di- 
mensions does not differ qualitatively from the behavior 
of g a (e — /i) in two dimensions. Some quantitative dif- 
ferences observed arise from the fact that at given N 
(the parameter which determines the amount of com- 
puter memory needed for the calculations) the linear size 
of a two-dimensional sample with a given density of sites 
is larger than that of a three-dimensional sample with the 
same density of sites and thereby, the finite size effects 
for three-dimensional samples with given N are more pro- 
nounced compared to those for the two-dimensional sam- 
ples with the same N. For example, the lower boundary 
of the region where g a (s — fi) can be described by the 
power law \e — fj,\ J shifts towards larger |e — fi\ > 0.4 
values (see inserts in Fig.|]). Remarkably, the exponent 
7 does not reach the value D — 1 predicted by the uni- 
versality hypothesis (|l|) even at A = (Fig.||). 

Unlike g a (s—fi) in the vicinity of the Coulomb gap, the 
density of ionized acceptors C a (|ll]) describes the state 
of the entire sample and therefore reaches the thermody- 
namic limit much faster than g a (s — fi). This allows us 
to obtain quite accurate results for C a from data on a 
relatively small amount of samples with TV = 500 only. 
Fig.^ shows the variations of C a with A both for two and 
three dimensions. In three dimensions almost all accep- 
tors become ionized (C Q ~ 1) rather soon while for two 
dimensions even for the largest A investigated around 10 
% of the acceptors remain neutral. So, one can say, that 
the three-dimensional MCDAM at A > 7 reduces already 
to the classical d— a model. It is known that the classical 
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TABLE I. The means fi and standard deviations A/j, of the 
Fermi energy calculated for the three-dimensional model ([]) 
with N = 1000 and various A. 



Secondly, as follows from the ground-state stability re- 



/' 



A [i 




2 
4 
6 
8 
10 



-0.017 
-1.016 
-2.0149 
-3.016 
-4.011 
-5.018 



0.100 
0.187 
0.287 
0.392 
0.502 
0.607 



d — a model exhibits in three dimensions the, so called, 
Coulomb fluctuational catastropheEj. For calculations 
on finite samples it implies that statistical fluctuations 
of /x(R-) grow dramatically with increasing A which is 
the case in our calculations (see Table |). Therefore, in 
order to reduce the statistical noise in three dimensions, 
the average of g a (e — jj) over a much larger (compared to 
D = 2) number of samples is needed. Note, that /x(R) in 
both two and three dimensions are scattered according to 
the Gaussian distribution with the mean p, obeying the 
relation (20). 



V. DISCUSSION 

The behavior of g a (e — /i) calculated within the re- 
gion of the Coulomb gap for the model @ is in strong 
contradiction to the universality hypothesis (Q) . Despite 
the fact that g a (s — fi) is indeed described by the power 
law \e — /i| 7 in a wide range of e inside the region of 
the Coulomb gap, the exponent 7 is considerably smaller 
than that predicted by the hypothesis (Q) both for the 
two- and three-dimensional cases. Moreover, the expo- 
nent 7 depends significantly on A and is different for the 
cases e > /i and e < \i. It is believed that information 
about g{e) might be directljj-abtained from tunneling aad. 
photoemission experimentsEa and recent experimentsE2l 
on boron-doped silicon crystals have shown that the den- 
sity of one-electron excitations at higher energies obeys a 
power-law with an exponent slightly less than 0.5 which 
is in good agreement with our results for D = 3 and 
A > 8. However, the non-metallic samples show around 
the Fermi energy a nearly quadratic Coulomb gap, so 
the question arises whether our results could be related 
to the intermediate asymptotic behavior observed? Here 
we want to make three remarks concerning this question: 

First, the power law g a (e — /x) ~ |e — /i| 7 is valid above 
a value So(N) below which the finite size effects take over 
(Figs -H and f§). It seems from our results, that eo(N) — > fi 
when N — > 00. In two dimensions we were able to obtain 
size- independent results down to £0 ~ 0.1, i. e. for ~ 90% 
of the whole Coulomb gap, the halfwidth of which is ~ 1 
in units of iso- 



lations 



the distance r,- 9 - between a neutral donor, with 



an energy, say, e\ G [—e,0] (e here is the halfwidth of a 
narrow band around /i = 0) and a charged donor with 
an energy G [0,e] should be not less that I. c., 
sites with energies e\ G [~£>0] cannot be inside a D- 
dimensional sphere of radius R sp = ^ and with the cen- 
ter in a site with the energy e° G [0, e]. Assuming that all 
such spheres do not intersect, the total volume occupied 
by the spheres is 



V sp = N X S{D) 



D 



g(e')de' 



(24) 



where S(D) is the volume of a D-dimensional sphere with 
the radius equal to unity. Since V sp cannot exceed the 



total volume V of a sample (V 
at the inequality 



N at n = 1) we arrive 



^ <- W 



which is valid for all e if 

D x 2 D 



(e) < 



\D-1 



S(D) 



(25) 



(26) 



le. 



The universality hypothesis (Q) then is a limit case of 
(p6|). The density of sites with energies e\ G [— e, 0} 
indeed decreases when e — > 0, so the assumption 
for the spheres with finite radii seems to be plausi 
However, simultaneously R sp — * 00 and consequently the 
plausibility of the assumption (pi| ) and thereby of the 
hypothesis (|l|) becomes questionable. 

And finally, the universality hypothesis (||) can be also 
obtained as the asymptotic behavior of a non-linear in- 
tegral equation for g(e) as e — > 0, the equation which, in 
turn, is heuristically obtained from the stability condition 
(Q). The derivation of this integral equation (given, for 
example, in Ref. |l7|) is based on the implicit assump- 
tion that the sites with charged donors are randomly 
distributed in space according to the Poisson statistics. 
However, it was unequivocalbz-jdemonstrated in computer 
studies of the Coulomb gaptiJ that charged donor sites 
with energies close to fi tend to form clusters (Ref. O, 
Fig. 6). 

We conclude that g a (e—fj,) in the region of the Coulomb 
gap in model (||) has a power law behavior for all energies 
down to [i and that the universality hypothesis of Efros 
(|l|) is questionable. Note, that our results are in contra- 
diction not only to the universality hypothesis (|l]), but 
to the inequality ( |26| ) as well. Up to now, all exponents 
found are in good agreement with this inequality. E. g. in 
Ref. [l2| specimens of 40 000 and 125 000 sites for two- and 
three-dimensional samples were investigated in the Efros' 
lattice modelllj and the power law g a (e — fi) ~ \e — /i| 7 
was found with 7 = 1.2 ± 0.1 and 7 = 2.6 ± 0.2 for two 
and three dimensions, respectively. The main conclusion 
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TABLE II. Some donor-acceptor pairs for which the difference between the donor and acceptor energy levels does not exceed 
10 meV. Eg, E v and E c are, respectively, the energy gap, the top of the valence band and the bottom of the conductivity band. 
If the solubilities of both donor and acceptor are known, the parameter Eo is calculated using the data for the less soluble of 
the pair. 



Donor 


Acceptor 


Solubility, cm 3 Eo, meV 
min max min max 


Ej, meV 


A, meV 








Si {E g = 1124 meV, x = 12) 






Fe 


Zn 


1.2 x 10 16 

2.3 x 10 16 


4 x 10 16 0.6 4 
8 x 10 16 


E c - 796 
E v + 320 


8 


Ni 


In 


10 18 
3 x 10 17 


10 13 12 < 1 
4 x 10 18 


E v + (160 -r 190) 
E v + 156.9 


3.1 - 33.1 








Ge {Eg = 740 meV, x = 15-9) 






S 


Ni 


no reliable data 
4.8 x 10 15 8 x 10 15 1.5 1.8 


E c - 296 
E c - 300 


4 








GaAs {E g = 1520 meV, x = 12-5) 






Ti 


Fe 


2 x 10 16 


3.1 


E c - 1000 
K + 520 






of ourpjjesults and those of Ref. |l2j is that Efros' lattice 
modellj can not be used to as a reliable approximation 
to the classical d — a model with Poisson impurity dis- 
tribution. Energy levels of donor (acceptor) impurities 
are usually close to the bottom (top) of the conduction 
(valence) band. Since in the most common semiconduc- 
tors the energy gap E g ~ 10 4 K and E ~ 20 K, A > 1 
and one may ask what physical relevance does the model 
(||) with a finite A < 10 have, except for being a pure 
academic exercise? However, in the case of deep impuri- 
ties the energy levels for some donor-acceptor pairs are 
extremely close to each other not excluding even the case 
A = OK Table [n] shows some donor-acceptor pairs with 
A < 10 in the most common semiconductors. The sol- 
ubilities of these impurities are rather low, thereby re- 
ducing the temperature at which the Coulomb gap with 
features described by the model (||) can be observed. For- 
tunately, these temperatures are high enough (~ 10 -j- 20 
K) for modern experimental techniques and hence exper- 
imental observation of the Coulomb gap in the semicon- 
ductors with deep impurities is possible to accomplish. 



VI. SUMMARY 

We have studied a model of impurities in semiconduc- 
tors with infinite-range Coulomb interactions between 
donors, between acceptors and between donors and ac- 
ceptors. A new parameter introduced in the model is 
the finite energy A of charge transfer between donors 
and acceptors, a parameter which enables processes of 
ionization of neutral impurities and of recombination 



of charged impurities. In the particular case of equal 
amounts of donor and acceptor impurities, we derived 
rigorous relations for the symmetry of the model with 
respect to exchange of donor and acceptor sites. We 
also extended the previously known algorithm to find the 
ground state including the stability relations with respect 
to ionization and recombination processes and performed 
computer studies of the model proposed at zero tempera- 
ture on a number of two- and three-dimensional samples 
with randomly distributed N donors and N acceptors. 
We explored the energy region around the Fermi energy 
/i where the Coulomb gap in the density of one-electron 
excitations g{e) is observed. The analysis of the calcu- 
lated histograms g{e) revealed that the behavior of g{e) 
obtained from the simulations on finite samples in the 
immediate neighborhood of /i is determined solely by the 
finite size effects. In the region where finite size effects 
become negligible g{e) is described by a power law with 
an exponent considerably depending on the parameter A 
and on the sign of e — [i. Our findings challenge the Efros 
universality hypothesis. Moreover, our results are in con- 
tradiction to the main inequality ( p6[ ) of which Efros' uni- 
versality hypothesis is a particular case. We have reex- 
amined the heuristic derivation of the Efros hypothesis 
and shown that some implicit assumptions which lead 
to universality are questionable. From the analysis of 
experimental data on admixtures in semiconductors we 
put forward possible experimental situations where one 
could observe the Coulomb gap with the features being 
the same as those of the model with a finite A. 
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